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Abstract 

A stochastic queueing network model with parameter-dependent 
service times and routing mechanism, and its related performance mea- 
sures are considered. An estimate of performance measure gradient is 
proposed, and rather general sufficient conditions for the estimate to 
be unbiased are given. A gradient estimation algorithm is also pre- 
sented, and its validity is briefly discussed. 

Key- Words: queueing networks, parameter-dependent routing, per- 
formance measure gradient, unbiased estimate. 

1 Introduction 

The evaluation of performance measure gradient presents one of the main 
issues of analysis of queueing network performance. Except in a few partic- 
ular models, there are generally no closed-form representations as functions 
of network parameters available for performance measures and their gradi- 
ents. In this situation, one normally applies the Monte Carlo approach to 
estimate gradient of network performance measures. 

In the last decade, infinitesimal perturbation analysis (IPA) [TJ [2J |BJ has 
received wide acceptance in queueing system performance evaluation as an 
efficient technique underlying the calculation of gradient estimates as well 
as the examination of their unbiasedness. Specifically, this technique was 
employed in [4] to calculate gradient estimates in closed networks with an 
ordinary probabilistic routing mechanism and general service time distri- 
butions. An extension of IPA, smoothed perturbation analysis, has been 
applied in [5] to the development of asymptotically unbiased gradient esti- 
mates in queueing networks with parameter-dependent routing. 

Another approach based on the analysis of algebraic representation of 
queueing system dynamics and their performance has been implemented 
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in O (3 [8]. This approach offers a convenient and unified way of analyt- 
ical study of gradient estimates, and it leads to computational procedures 
closely similar to those of IPA. In this paper, based on this approach, a 
rather general queueing network model with parameter-dependent service 
times and routing mechanism is presented. For the performance measures 
which one normally chooses in analysis of network performance, we propose 
a gradient estimate, and give sufficient conditions for the estimate to be un- 
biased. These conditions are rather general and normally met in analysis of 
queueing network performance. Finally, an algorithm of estimating gradient 
of a particular performance measure is presented, and its validity is briefly 
discussed. 

2 The Underlying Network Model 

We consider a generalized model of a queueing network consisting of N 
nodes, with customers of a single class. As is customary in queueing network 
models, customers are assumed to circulate through the network to receive 
service at appropriate nodes. We do not restrict ourselves to a particular 
type of nodes, it is suggested that any node may have a single server as well 
as several servers operating either in parallel or in tandem. 

Furthermore, there is a buffer with infinite capacity in each node, in 
which customers are placed at their arrival to wait for service if it cannot 
be initiated immediately. We assume the queue discipline underlying the 
operation of any node to be first-come, first-served. Upon his service com- 
pletion at one node, each customer goes to another node chosen according 
to some routing procedure described below. We suppose that the transition 
of any customer between nodes requires no time, and he therefore arrives 
immediately into the next node. Finally, we assume that the network starts 
operating at time zero; at the initial time, the server at any node n is free, 
whereas its buffer contains K n customers, < K n <oo, n = 1, . . . , N . 

We now turn to the formal description of the network dynamics from 
an algebraic viewpoint, and then introduce randomness into the network 
model. 

2.1 Algebraic Description of Node Dynamics 

In a general sense, each node can be regarded as a processor which produces 
an output sequence of departure times of customers from another two, an 
input and control sequences formed respectively by the arrival times and the 
service times of customers. Let us denote for every node n, 1 < n < N , the 
A;th arrival epoch to the node by A\ , and the fcth departure epoch from the 
node by D k . Notice, because the transition of customers from one node to 
another is immediate, each A 1 ^ coincides with some Dj with the exception 
of A k n = for all k < K n . Finally, we denote the service time associated 



2 



with the kth. service initiation in node n, by r*. The set of all service times 
T = {r k \n = 1, . . . , N; k = 1, 2, . . .} is assumed to be given. 

The usual way to represent the operation of a node is based on recursive 
equations describing evolution of D k as a state variable [9J \W[ [11] . Note 
that these recursive equations are often rather difficult to resolve. Below 
are given two equations which describe dynamics of nodes operating as the 
G/G/l and G/G/2 queueing systems. Other examples may be found in 

u MM, mj. 

2.1.1 The G/G/l queue. 

Suppose first that node n is represented as the G/G/l queue. Its associated 
recursive equation may be written as [9] 

D> = iA*VD*T 1 )+T* i 

where V denotes the maximum operator, and D k = for all k < 0. It is 
easy to see that the solution of the equation in terms of arrival and service 
times has the form 

k I k 

D k n = v [<+Y, T i 

i=l y j=i 

2.1.2 The G/G/2 queue. 

The equation which describes the dynamics of a node operating as the 
G/G/2 queue may be considered as rather difficult to handle. For node 
n, it is written as [12] 

k 

D k n = M ((4 V Z?r 2 ) + <) A [(A k n +1 V Dt 1 ) + r n fc+1 ) , 

i=l 

where A stands for the minimum operator. Although there are no closed- 
form solutions of the equation, known to the author, it is clear that it exists. 



2.2 Routing Mechanism and Interaction of Nodes 

The routing mechanism inherent in the network is defined by the sequences 
R n = {Pn,Pni • • •} given for each node n, where p\ represents the next 
node to be visited by the customer who is the &th to depart from node n, 
p k n G {1, . . . , N}, k = 1, 2, . . . . The matrix 
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is referred to as the routing table of the network. 

In order to describe the dynamics of the network completely, it remains 
to define formally interactions between nodes. In fact, a relationship between 
arrival and departure times of distinct nodes is to be established. To this 
end, for each node n, let us introduce the set 

D n = {Dj\4 = n; i = 1, . . . , N; j = 1, 2, . . .} (1) 

which is constituted by the departure times of the customers who have to go 
to node n. Furthermore, we denote by A\ the arrival time of the customer 
which is the /cth to arrive into node n after his service at any node of the 
network. In other words, the symbol A\ differs from A k in that it refers 
only to the customers really arriving into node n, and does not to those 
occurring in this node at the initial time. 
It has been shown in 18] that it holds 

A k n = f\ piV-V^), (2) 

{Di,...,D fc }C.D„ 

where minimum is taken over all A: -subsets of the set D n . The times A\ 
and A\ are related by the equality 

A k = { °' , if k ~ Kn (3) 
n \ An~ Kn , otherwise ' ^ ' 

Clearly, if deterministic routing with an integer matrix R as the routing 
table is adopted in the model, each set D n , n = 1, . . . , N, is determined 
uniquely from ([1]), and then a straightforward algebraic representation of A^ 
may be obtained using (J2IEI) • In this case, starting from the above represen- 
tations of node dynamics, one may eventually arrive at algebraic expressions 
for any arrival time A^ and departure time D k , which are written in terms 
of service times r £ T , and involve only the operations of maximum, mini- 
mum, and addition. 

2.3 Representation of Network Performance 

One of the features of the formal network model described above is that it 
offers the potential for representing network performance criteria in a rather 
simple and convenient way. Suppose that we observe the network until the 
i^th service completion at node n, 1 < n < N . As performance criteria 
for node n in the observation period, one normally chooses the following 
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average quantities El El El [JT] : 



System time of one customer, 


OA _ 
D n — 






Waiting time of one customer, 


W K = 




-A k n -r k )/K, 


Throughput rate of the node, 


rpK _ 


K/D*, 




Utilization of the server, 


U K = 




K 
n ) 


Number of customers, 


J K ~ 




A n )/D n , 


Queue length at the node, 


Q« = 







It is easy to see that with the routing mechanism determined by an integer 
matrix, all these criteria may be represented only in terms of service times 
in closed form. 

2.4 Stochastic Aspect and Performance Evaluation 

Let us suppose that for all n = 1, . . . , N , and k = 1,2, . . . , the service times 
are defined as random variables r k = T k (6,u), where 6 G 6 C K is a deci- 
sion parameter, w is a random vector which represents the random effects 
involved in network behaviour. First we assume the routing table R = R to 
be an integer matrix. Since deterministic routing leads to algebraic expres- 
sions in terms of the random variables r G T for the performance criteria 
introduced above, one can conclude that these criteria also present random 
variables. 

Let F = F(6, uj) be a random performance criterion of the network. As 
is customary, we define the performance measure associated with F as the 
expected value 

F(0)=E w [F(9,lj)]. (4) 

Although we may express F in closed form, it is often very difficult or impos- 
sible to obtain analytically the performance measure F . In this situation, 
one generally applies a simulation technique which allows of obtaining values 
of F(9,u) , and then estimate the network performance by using the Monte 
Carlo approach. 

We now turn to the networks with parameter-dependent probabilistic 
routing. We assume p k = p„(0,u;) to be a discrete random variable ranging 
over the set {1, . . . , N} . The routing mechanism of the network is now 
defined by the random matrix R = R(9,uj) with particular routing tables 
as its values. We denote the set of all possible routing tables R by 1Z. 

Obviously, the expression of A k n defined by fl2ED 

may change from one 

shape into another, depending on particular routing tables. To take this 
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into account, we now define the random performance criteria in j3D as EE] 

F(e,u, R{e,u)) = Y, ifflfM^^t^)- ( 5 ) 

Ren 

where 1{h(6»,o;)=_r} is the indicator function of the event {R(9,u) = R} , 
and Fr(6,uj) = F(9,uj,R) is the performance criterion evaluated under the 
condition that the network operates according to the deterministic routing 
mechanism defined by the routing table R € 1Z. 



3 Performance Measure Gradient Estimation 

Since there are generally no explicit representations as functions of system 
parameters 9 , available for the performance measure F , one may evaluate 
its gradient dF(9)/d9 by no way other than through the use of either finite 
difference estimates [HI E] or the estimate 

1 M 

g(6,uj 1 ,...,uJM) = Jj^2-QQ F ( 9 ^i)^ ( 6 ) 

8=1 

where , i = 1, . . . , M , are independent realization of oj , provided that the 
derivative dF(9,u)/d9 exists. 

Very efficient procedures of obtaining gradient estimates may be designed 
using the IPA technique [H El [3]. Such a procedure can yield the exact 
values of the derivative dF(9, u)/d9 by performing only one simulation run. 
Furthermore, in the case of a vector parameter 9 € R d , the IPA procedures 
provide all partial derivatives dF(9,uj)/d9i , i = l,...,d, simultaneously, 
and take an additional computational cost which is normally very small 
compared with that required for the simulation run alone. Finally, it can 
be easily shown [131 [3] that if the IPA estimate of the derivative in ([6]) is 
unbiased, the mean square error of g has the order which is significantly 
less than those of any finite difference estimates based on the same number 
of simulation runs. 

A sufficient condition for the estimate (|6j) to be unbiased at some 9 E 
requires [13l El [3] 



^E[F(9,u)]=E 



(7) 



A usual way of examining the interchange in ([7]) involves the application of 
the Lebesgue dominated convergence theorem 



Theorem 1. Let (U,F,P) be a probability space, & C \R d , and F : QxQ — > 
R be a T -measurable function for any 9 € O and such that the following 
conditions hold: 
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(i) for every 9 £ Q, there exists dF(9,uj)/d9 at 9 w.p. 1, 

(ii) for all 61,62 £ O, there is a random variable X(uj) with EX < oo and 

such that 

\F(9 1 ,u)-F(9 2 ,u)\<X(u))\\9 1 -e 2 \\ w.p. 1. (8) 

Then equation (|7|) holds on . 

In [6l [8] the approach based on the implementation of Theorem Q] has 
been applied to analyze estimates of performance gradient in the networks 
models with the parameter-independent probabilistic routing mechanism de- 
termined by a random routing table R = R(u). Specifically, starting from 
the representations of network dynamics, discussed in Section 2, it has been 
shown that 

(i) if each service time r G T satisfies the conditions of Theorem [IJ and for 

every 9 £ O, all r £ T present continuous and independent random 
variables, then the average total time and waiting time satisfy 
the conditions of Theorem [IJ 

(ii) if in addition to previous assumptions, there exist random variables 
fx, v > such that for all r £ T it holds v < |r| < \i w.p. 1 for all 
9 £ @, and the condition E[fiX/u 2 ] < oo is fulfilled, where A is the 
random variable providing r with ((8|), then the average throughput 
rate T„ , utilization U*~ , number of customers , and queue length 
Qn satisfy the conditions of Theorem [TJ 

Note that the above conditions do not involve independence at each 9 be- 
tween the random variables r(9, u) and p[oS) in the probabilistic sense. 

4 An Unbiased Gradient Estimate for Networks 
with Parameter-Dependent Routing 

We start the section with an example which exhibits difficulties arising in 
gradient estimation when there is a parameter dependence involved in the 
routing mechanism of the network, and then present our main result offering 
an unbiased estimate of performance measure gradient in networks with 
parameter-dependent routing. 

4.1 Preliminary Analysis 

Suppose that there is a parameter dependence of the routing mechanism in 
the network, that is R = R{9,uj). In this case, the random performance 
criterion F generally violates condition (JH]). As an illustration, one can 
consider the following example. 
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Let 6 = [0, 1], Q\ = f2 2 = [0, 1], and (f2, F, P) be a probability space, 
where F is the a -field of Borel sets of f2 = Sl\ x f2 2 , -P is the Lebesgue 
measure on Q. Denote oj = (oji,oj 2 ) t , and define the function 



F[e,u, P [e,u)) 

where 



n(e,uj), ifp(0,w) = l 
t 2 {0,uj), if p(0,w) = 2 ' 



and 

r i, ifw 2 <0 
^' w) = \ 2, if. 2 >^ • 

We may treat t\ and t 2 as the service time of a customer respectively at 
node 1 and 2. The function F is then assumed to be the service time of 
the customer which may arrive into either node 1 or 2, according to one of 
the two possible values of p. 

Clearly, t\ and r 2 satisfy the conditions of Theorem [IJ whereas the 
function F now represented as 



F(6,u>) 



+ U1X + 1, if w 2 < # 
+ uj\, if uj 2 > 9 ' 



is differentiable w.p. 1 at any 8 € G, and dF(6,u>)/d8 = 1 w.p. 1. However, 
for any 0\ , 6 2 € G such that 9\ > uj 2 and 6 2 < oo 2 , it holds 

\F(e u co) -F(e 2 ,u) | > i, 

and therefore, condition ([8|) is violated. 

On the other hand, it is easy to verify that 

E[F(0, oj)] = 26 + |, &E[F{0, oj)) = 2 

^F(0,a;) = l w.p. 1, £7[$F(0,a;)] = 1. 

In other words, equation (|7|) proves to be not valid, and we finally conclude 
that estimate (0) will be biased. 

4.2 The Main Result 

To suppress the bias in estimates of the gradient 

^F{6) = ^E[F(6,u,R(d,cj))], (9) 
let us replace © by the estimate 

1 - 

g(9,u) 1 ,...,u M ) = —Y, G (^)< ( 10 ) 
where G(9,u>) will be defined in the next theorem. 
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Theorem 2. Suppose that a random performance criterion F is represented 
in form ([5]) , and for each R € 1Z the following conditions hold: 

(i) Fr satisfies the conditions of Theorem^ 

(ii) for any 9 6 0, the random variable Fr(0,ui) and the random matrix 

R(9,lu) are independent; 

(iii) for any S £ 0, the function &(0,R) = Pr{R(9,co) = R} is continu- 
ously differ entiable at 9, and &(9,R) > 0. 



Then for any 9q £ Q , estimate (fTU|) with 

+ F(e ,u,R(e ,uj))y(9o,R{e ,u;)), 



G(9 ,u) = ^F(9,u,R(9 ,uj)) 



e=e 

where *$>(9,R) = In &(9, R) , is unbiased. 

Proof. Clearly, it is sufficient to show that the equation 

^E[F(9,u,R(9,cj))]=E[G(9,gj)) 

holds for any 9 £ 0. 

To verify this equation, let us first represent F in form ([5|), and consider 
its expected value 

E[F(0,u, R(6,w))} =Y, E lHR{e,u)=R}F{9,u,R)] . 
Ren 

Since E[liRfg u \ = m) = Pt{R(9,lj) = R} = &(9,R), it follows from condi- 
tion (ii) of the theorem that 

E[F(0 t u t R(9,u))] 

= E[F(9, u, R)]Pt{R(0, lu) = R}= ^ E[F(8, u, R)}<&(9, R). 
Reii Ren 
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For any 9q £ , under conditions (i) and (iii) , we successively get 

^E[F(e,u,R(e,oj))] 



8=9o 



E sW^i 



H0 ,R) 



e=e 







+ E[F(9 ,lo,R)} —$(9,R) 



E \ E 



Ren 



d_ 
dd 



F(9,u>,R)] 



9=9 J 



e=e 



+ E[F(e ,u,R)] 



${0 ,R) 



E * 

Ren \ 



d_ 
89 



E E 

Ren 



d_ 



F(d,u,R) 



R) 



= 00 



+ E[F(6 ,L,,R)]*(e ,R) )$(9 ,R) 



+ F(9 q ,lo,R)^(0 o ,R) 



Pt{R(9 ,u) = R} 



E 



d 



F(9,u,R(9 ,lu)) 



+ F(9 ,u,R(9 ,u))^(9 ,R(9 ,cu)) 



8=9o 



09 

= E[G(9 ,uj)). □ 

It is not difficult to obtain the conditions for estimate (|10p to be unbi- 
ased for gradient of particular performance measures. They can be stated 
by combining the conditions in Section 3, related to particular performance 
criteria in networks with parameter-independent routing, with those of The- 
orem [2j Note that these conditions are rather general, and normally met in 
analysis of queueing networks. 

Let us now return to the example presented in the previous subsection. 
First, we have 

$(9,l)=Pv{p(9,uj) = l} = 9, 
<S>(9,2) = Pr{p(9,cu) = 2} = l-0, 

and then 



d 



•(«,*) = 3 Mi 



1 



1 



In this case, the function G is defined as 

1 + «±«i±l i io j 2 <0 



G(9,uj) 
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for any 9 £ (0, 1) . Finally, evaluation of its expected value gives 
E[G(9,lo)] = ^E[F(8,lj,p(9,lj))] = 2. 

5 Application to Network Simulation 

Consider a network with N single-server nodes, and assume that the Kth 
service completion at a fixed node n, 1 < n < N , comes with probability 
one after a finite number of service completions in the network. In this case, 
to observe evolution of the network until the i^th completion at the node, 
it will suffice to take into consideration only finite routes defined by a right 
truncated routing table with integer (N x L)-matrices 



R = 


( 


r\ 


r\ . 
r\ . 


. r\ 
■ r\ 


\ 




\ 


r N 


r N ■ 


■ r N 


J 



as its values, with some L > K . 

Furthermore, we assume, as is customary in network simulation, that 
for each 6 £ 0, the random variables p!^(9,uj) are independent for all 
n = 1, . . . , N , and k = 1, . . . , L. With this condition and the notation 
9^(6*, r) = ~Pr{p k n {6, u) = r}, we may represent the function $ introduced 
in Theorem [21 as 

N L 
n=l k=l 

and then get the function ^/ in the form 

8 N L 8 

9(0, R) = gE^ HO, R) = E E OS ln 

n=l k=l 

Suppose now that we have to estimate the gradient of a performance 
measure, say U%(6) = E[U^(9,oj)}, the expected value of the average 
utilization of the server at node n. It results from Theorem [2] that, as an 
unbiased estimate, the function 

G(6,u>) = ^F(e,u,R) + F(e,u,R)9(9,R), (11) 
may be applied with R = R(0,u), and 

K 

F(9,u 1 ,R) = J2'4(0,oj)/D^(e,cj,R). 
k=i 



11 



It is not difficult to construct the next algorithm which produces the 
value of G(8,uS) for fixed 8 £ C R, and u € f2, provided that there is 
a network simulation procedure into which the algorithm may be incorpo- 
rated. It actually combines an IPA algorithm [1] for obtaining the derivative 
8F(8,uj,R)/88 with additional computations according to (jlip . 

Algorithm 5.1. 

Initialization: 

for i = 1, . . . , N do gi i — ; 

s,t,t' i — 0; 
R < — R(0,lo); 



Upon the kth service completion at node i, perform the instructions: 

if i = n then t < — t + t^(9,lj); 

if k = K then d < — D% (0, u, R) ; 
stop; 

r < — r 1 ; 

«<— a + ,|]n^(0,r); 

if the server of node r is free then g r < — . 

On completion of the algorithm, it remains to compute (t'd — tg n )/d 2 + 
ts/d as the value of G. 

Note, in conclusion, that estimate ([6]) with the function G evaluated 
using the algorithm will not be unbiased in general. For the estimate to 
be unbiased, the function ^(9,R) in (jlip must be calculated as the sum 
with the same number of summands d In ip^(9 , r) / 88 for any of simulation 
runs. However, during the simulation runs with distinct realizations of oo, 
there may be different numbers of service completions encountered at nodes 
i 7^ n, and then considered in evaluation of . This normally involves an 
insignificant error in estimating the gradient, which becomes inessential as 
K increases. 
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